3 





Discretization 



Discretization takes the fundamental idea of calculus 
and pushes it to the opposite extreme from what cal- 
culus uses. Calculus was invented to analyze chang- 
ing processes such as orbits of planets or, as a one- 
dimensional illustration, how far a ball drops in time 
t. The usual computation 

distance = velocity x time 

fails because the velocity changes (it increases linearly with time). However 
- and this next step is the fundamental idea of calculus - over a short time 
interval, its velocity is almost constant and the usual distance formula works 
for each short interval. Each short distance is the area of one rectangle, and 
the total distance fallen is approximately the combined area of the rectangles. 
To eliminate the error, calculus uses the extreme case of infinite rectangles, 
ever thinner (shorter intervals) until each shrinks to zero width. Then the 
approximation of constant speed becomes exact. Discretization uses the op- 
posite extreme: one maybe two fat rectangles. This limitation means the 
error may not be zero, but it drastically simplifies any computations. 




3.1 Exponential decay 

The first example is this integral: 1 

/>oo 

/ e~* dt. 
Jo 

Since the derivative of e* is e*, the indefinite integral o 
is easy to find exactly, and the limits make the compu- 
tation even simpler. In an example where the exact answer is known, we can 
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adjust the free parameters in the method of discretization until the method 
produces accurate values. So, replace the complicated, continuous, smooth 
exponential decay e~* by a rectangle, and do the integral by finding the area 
of the rectangle. 



With one rectangle, the approximate function re- 
mains constant until it abruptly falls to and remains 
zero. Finding the area of the rectangle requires choos- 
ing its height and width. A natural height is the max- 




imum of e _t , which is 1. A natural width is the time Q _| 

interval until f(t) = e _t changes significantly. A sim- 1 
ilar idea appeared in Section 1.4 to approximate a 

derivative df jdx. Its numerator df was estimate as a typical value of f(x). 
Its denominator dx became the x interval over which f(x) changes signifi- 
cantly. For an exponential, a natural definition for significant change is to 
changes by a factor of e. When f(t) = e~* , this change happens when t goes 
from t to t + 1. So the approximating rectangle, whose height we've chosen 
to be 1, also has unit width. It is a unit square with unit area. And this 
rectangle exactly estimates the integral since 

poo 

/ e~ t dt=l. 
Jo 



3.2 Circuit with exponential decay 



In Chapter 1 on dimensions, I insisted that declaring quantities prematurely 
dimensionless ties one hand behind your back. In the previous example I 
committed that sin by making the exponent be —t. Since an exponent is 
dimensionless, my choice made t dimensionless as well. 

A more natural interpretation of t is as a time. So here is a similar 
example but where t has dimensions, which are useful in making and 
checking the approximations. Let's first investigate the initial condi- 
tions, just before the switch closes. No current is flowing since the 
circuit is not yet a closed loop. Furthermore, because the circuit has ( 
been waiting forever, the capacitor has had completely discharged. 
So capacitor has no charge on it. The charge determines the voltage 
across the capacitor by 



Q = CV C 
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where Q is the charge on the capacitor, C is its capacitance, and Vc is the 
resulting voltage. [See the classics on circuits [2] and electromagnetism [3] for 
more on capacitors.] So just before the switch closes, the capacitor has zero 
voltage on it (Vc = 0). 

At time t = 0, 1 close the switch, which connects the resistor and capacitor 
to the source voltage V (which is constant). Since Vc starts at zero, the 
voltage drops in the resistor is the whole source voltage V: 

V R = V (initially), 

where Vr is the voltage across the resistor. This voltage drop is caused 
by a current / flowing through the resistor (which then flows through the 
capacitor). Ohm's law says that Vr = IR. Initially Vr = V so the initial 
current is lo = V/R. This current charges the capacitor and increases Vc- As 
Vc increases, Vr decreases - which decreases the current /, which decreases 
how fast Vc increases, which . . . Finding the current is a problem for calculus, 
in particular for a differential equation. 

Instead, let's guess the current using dimensions, 

, , i r , I(t->Q) 

extreme cases, and the new technique of discretiza- Io — i 1 

tion. First apply extreme cases. At the t = ex- 
treme, the current is I — V/R. At the t = oo ex- 
treme, no current flows: The capacitor accumulates 

enough charge so that Vc = V, whereupon no voltage o — j — — - — — ;: 

drops across the resistor. From Ohm's law again, a 
zero voltage drop is possible only if no current flows. 

Between those extremes, we guess / using discretization. Pretend that / 
stays at its t — value of Iq for a time r, then drops to its t = oo value of 
1 = 0. So t is the time for the current to change significantly. To determine 
r, use dimensions. It can depend only on R and C. [It could depend on 
V, but because the system is linear, the time constants do not depend on 
amplitude.] The only way to combine R and C into a time is the product 
RC. A reasonable guess for r is therefore r = RC. In this picture, the 
discretized current stays at V/R until t = r, then falls to and remains zero 
forever. 
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discretized / 



No physical current changes so abruptly. To guess 
the true current, use discretization in reverse. The 
exponential decay of Section 3.1 produced the same 
rectangular shape after discretizing. So perhaps the 
true current here is also an exponential decay. In 
the other example, the function was e _t , and the 
changeover from early- to late-time behavior happened 

at t = 1 (in that example, t had no dimensions). By t = 1, the exponential 
decay e~* had changed significantly (by a factor of e). For this circuit, the 
corresponding changeover time is r. To change by a factor of e in time r, 
the current should contain e~*/ T . The initial current is / — Iq, so the current 
should be 




t/r 



Having a solution, even a guess, turns the hard work of solving differential 
equations into the easier work of verifying a solution. 

To test the guess for /, I derive the differential equation for the current. 
The source voltage V drops only in the resistor and capacitor, so their voltage 
drops must add to V: 



V = V R + V c . 

The capacitor voltage is Vc = Q/C. The resistor voltage is Vr = IR, so 

Q 



V = IR 



C 



It seems that there are too many variables: V and C are constants, but / and 
Q are unknown. Fortunately current / and charge Q are connected: charge 
is the time integral of current and / = dQ/dt. Differentiating each term with 
respect to time simplifies the equation: 








R 



dl 

dt 



I 

C' 



Move the R to be near its companion C (divide by R): 

dl I dl I 

di + ~RC ~ di + ' r' 



= 



Dimensions, extreme cases, and reverse discretization produced this current: 
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Amazing! It satisfies the differential equation: 



dt 



{'"-"') 



= 



because the time derivative brings down a factor of — 1/r, making the first 
and second terms equal except for a minus sign. 



3.3 Population 



discrctized distribution 



10 6 /year 




Age (years) 



Not all problems are exponen- 
tial decays. In the next exam- 
ple, the true functions are un- 
known and exact answers are 
not available. The problem is 
to estimate the number of ba- 
bies in the United States. To 
specify the problem, define ba- 
bies as children less than two 
years old. One estimate comes 
from census data, which is ac- 
curate within the limits of sta- 
tistical sampling. You integrate the population curve over the range t = 
to 2 years. But that method relies on the massive statistical efforts of the 
US census bureau and would not work on a desert island. If only the pop- 
ulation were constant (didn't depend on age), then the integrals are easy! 
The desert-island, back-of-thc-cnvclopc method is to replace the complicated 
population curve by a single rectangle. 

How high is the rectangle and how wide is it? The width t, which is 
a time, has a reasonable estimate as the average life expectancy. So t ~ 
70 years. How high is the rectangle? The height does not have such an 
obvious direct answer as the width. In the exponential-decay examples, the 
height was the the initial value, from which we found the area. Here, the 
procedure reverses. You know the area - the population of the United States 
- from which you find the height. So, with the area being 3 • 10 8 , the height 
is 



height 



area 3 • 10 8 
width 75 years 
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since the width is the life expectancy, for which we used 70 years. How did 
it become 75 years? The answer is by a useful fudge. The new number 75 
divides into 3 (or 300) more easily than 70 does. So change the life expectancy 
to ease the mental calculations. The inaccuracies caused by that fudge are 
no worse than in replacing the complex population curve by a rectangle. So 

height — 4 • 10 6 year -1 . 

Integrating a rectangle of that height over the infancy duration of 2 years 
gives 

^babies ~ 4 • 10 6 years" 1 x 2 years = 8 ■ 10 6 . 
height infancy 

Thus roughly 8 million babies live in the United States. From this figure, 
you can estimate the landfill volume used each year by disposable diapers 
(nappies) . 



3.4 Full width at half maximum 



The Gaussian integral 




-1 



has appeared in several examples, and you've seen the trick (in 
Section 2.2) of squaring it to show that its value is \pK. The exponential 
in the integrand is a difficult, continuous function. Except over the infinite 
integration range, the integral has no closed form, which is why statistics 
tables enumerate the related error function numerically. I introduce that 
evidence to show you how difficult the integral is without infinite limits, and 
it is not easy even with infinite limits. 

Pretend therefore that you forget the trick. You can ap- 
proximate the integral using discretization by replacing the in- 
tegrand with a rectangle. How high and how wide should the 
rectangle be? The recipe is to take the height as the maximum 
height of the function and the width as the distance until the function falls 
significantly. In the exponential-decay examples, significant meant changing 
by a factor of e. The maximum of e~ x is at x = when e~ x = 1, so the 
approximating rectangle has unit height. It falls to 1/e when x = ±1, so the 
approximating rectangle has width 2 and therefore area 2. This estimate is 




2008-03-06 13:24:47 / rev ebd336097912+ 



Cite as: Sanjoy Mahajan, course materials for 18.098 / 6.099 Street- Fighting Mathematics, IAP 2008. 
MIT OpenCourseWare (http://ocw.mit.edu/), Massachusetts Institute of Technology. 
Downloaded on [DD Month YYYY]. 



3.4 Full width at half maximum 



37 



half decent. The true value is y/n = 1.77. . ., so the error is about 13%: a 
reasonable trade for one line of work. 

Another recipe, also worth knowing because it is sometimes more accu- 
rate, arose in the bygone days of spectroscopy. Spectroscopes measure the 
wavelengths (or frequencies) where a molecule absorbed radiation and the 
corresponding absorption strengths. These data provided an early probe into 
the structure of atoms and molecules, and was essential to the development 
of quantum theory [4]. An analogous investigation occurs in today's particle 
accelerators - colloquially, atom smashers - such as SLAC in California and 
CERN and in Geneva: particles, perhaps protons and neutrons, collide at high 
energies, showering fragments that carry information about the structure of 
the original particles. Or, to understand how a finely engineered wristwatch 
works, hammer it and see what the flying shards and springs reveal. 

The spectroscope was a milder tool. A chart recorder plotted the absorp- 
tion as the spectroscope swept through the wavelength range. The area of 
the peaks was an important datum, and whole books like [5] are filled with 
these measurements. Over half a century before digital chart recorders and 
computerized numerical integration , how did one compute these areas? The 
recipe was the FWHM. 



FWHM = full width at half maximum 



Unpack the acronym in slow motion: 

1. M. Find the maximum value (the peak value). 

2. HM. Find one-half of the maximum value, which is the half maximum. 

3. FWHM. Find the two wavelengths - above and below the peak - where 
the function has fallen to one-half of the maximum value. The full width 
is the difference between the above and below wavelengths. 

The FWHM approximation recipe replaces the peak by a rectangle with 
height equal to the peak height and width equal to the the width estimated 
by the preceding three-step procedure. 

Try this recipe on the Gaussian integral and compare the 
estimate with the estimate from the old recipe of finding where 
the function changed by a factor of e. The Gaussian has max- 
imum height 1 at jc = 0. The half maximum is 1/2, which 
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happens when x = ±\An 2. The full width is then 2\/ln 2, and 

the area of the rectangle - which estimates the original integral - is 2\An 2. 

Here, side by side, are the estimate and the exact integral: 

^ = 1.7724... (exact), 
2Vrn2 = 1.6651... (estimate). 

The FWHM estimate is accurate to 6%, twice as accurate as the previous 
recipe. It's far better than one has a right to expect for doing only two lines 
of algebra! 



3.5 Stirling's formula 

The FWHM result accurately estimates one of the most useful quantities in 
applied mathematics: 

nl = n x (n — 1) x (n — 2) x • • • x 2 x 1. 

We meet this quantity again as a picture proof in Section 4.6. Here we 
estimate nl by discretizing an integral representing nl: 

poo 

/ fe - ' dt = nl 

You may not yet know that this integral is nl; you can show it either with 
integration by parts or see ?? on generalization to learn differentiation under 
the integral sign. For now accept the integral representation on faith, with a 
promise to redeem the trust in a later chapter. 

To approximate the integral, imagine what the inte- 
grand i"e~* looks like. It is the product of the increasing 
function t n and the decreasing function e - *. Such a prod- 
uct usually peaks. A familiar example of this principle is 
the product of the increasing function x and the decreas- 
ing function 1 — x. over the range x £ [0,1] where both 
functions are positive. The product rises from and then falls back to zero, 
with a peak at x = 1/2. 

You can check that the product £™e -t has a peak by looking at its behavior 
in two extreme cases: in the short run t — > and in the long run t — > oo. 
When t — ► 0, the exponential is 1, but the polynomial factor t n wipes it 
out by multiplying by zero. When t — > oo, the polynomial factor t n pushes 
the product to infinity while the exponential factor e - ' pushes it to zero. 
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An exponential beats any polynomial. To see why and avoid the negative 
exponent —t muddying this issue, compare instead e* with t n as t — > oo. The 
Taylor series for e* contains all powers of t, so it is like an infinite-degree 
polynomial. So e*/i" goes to infinity once t gets large enough. Similarly, 
its reciprocal t n e~ l goes to zero as t — > oo. Being zero at also f = 0, the 
product is zero at both extremes and positive elsewhere. Therefore it peaks 
in between. Maybe it has more than one peak, but it should have at least 
one peak. Furthermore, as n increases, the t n polynomial factor strengthens, 
so the e~* requires a larger t to beat down the t n . Therefore, as n increases 
the peak moves right. 

With fe - ' having a peak, the FWHM recipe 
can approximate its area. The recipe requires find- 
ing the height (the maximum of the function) and 
the width (the FWHM) of the approximating rec- 
tangle. To find these parameters, slurp the t n into 
the exponent: 

^Jig-i gTllntg— t ^nlnt—t 




The exponent f(t) = n\nt — t is interesting. As t — > 0, the hit takes f(t) to 
— oo. As t — > oo, the —t takes f(t) again to — oo. Between these limits, it 
peaks. To find the maximum, set f'(t) = 0: 

/'(t) = ^-l = 0, 

or ipoak = n. As we predicted, the peak moves right as n increases. The 
height of the peak is one item needed to estimate the rectangle's area. At the 
peak, f(t) is f(n) = nlnn — n, so the original integrand, which is e^ l \ is 



o/(^pcak) 



J(n) 



nlnn—n 



n" 

e" 



(=)"• 



To find the width, look closely at how f(t) behaves near the peak t — n by 
writing it as a Taylor series around the peak: 

f(t) - f(n) + f'(n)(t -n) + \f"{n){t - nf + ■ ■ ■ . 

The first derivative is zero because the expansion point, t = n, is a maximum 
and there f'(n) = 0. So the second term in the Taylor series vanishes. To 
evaluate the third term, compute the second derivative of / at t = n: 
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So 



f(t) = nlnn 




n) 2 + 



f(n) 



f"{n) 



The first term gives the height of the peak that we already computed. The 
second term says how the height falls as t moves away from n. The result is 
an approximation for the integrand: 



„/(*) 



-(;)' 



-(t-n) 2 /2n 



v'8nln2- 



The first factor is a constant, the peak height. The 
second factor is the familiar Gaussian. This one is 
centered at t = n and contains l/2n in the expo- 
nent but otherwise it's the usual Gaussian with a 
quadratic exponent. It falls by a factor of 2 when 
(f — n) 2 /2n = In 2, which is when 

t±=n± V2nln2. 

The FWHM is t+ - t_, which is V8nln2. The estimated area under e /(t) is 
then 



[-) xV8nln2. 



As an estimate for nl, each piece is correct except for the constant factor. The 
more accurate answer has \/2tt instead of \/81n2. However, 2tt is roughly 
8 In 2 so the approximate is, like the estimate the vanilla Gaussian integral 
(coincidence?), accurate to 6%. 



3.6 Pendulum period 

The period of a pendulum is by now a familiar topic in this book. Its differen- 
tial equation becomes tractable with a bit of discretization. The differential 
equation that describes pendulum motion is 



d 2 g 



= 



This nonlinear equation has no solution in terms of the usual functions - 
to put it more precisely, in terms of elementary functions. But you can 
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understand a lot about how it behaves by discretizing. If only the equation 
were 



cP0 



= 0. 



This equation is linear, and therefore possible to solve without too much 
misery - I hesitate to say that any differential equation is 'easy' - and its 
solution is an oscillation with angular frequency to = \J g/l: 



9{t) = 6> cos (jtj 



Its period is 2iT^g/l, which is independent of amplitude 9q. 

The complexity of the unapproximated pendulum equation 
arises because it has the torque-producing factor sin 9 instead 
of its approximation 9. The two functions match perfectly as 
9 — > 0. But as 9 grows - which happens with large amplitudes 
- the equivalence becomes less accurate. One way to compare 
them is to look at their ratio (sin 9) /9. As expected, when 9 = 
0, the ratio is 1. As 9 grows, the ratio falls, making the simple- 
harmonic approximation less accurate. We can discretize to 
find a more accurate approximation than the usual simple- 
harmonic one, yet still produce a linear differential equation. The upcoming 
figures illustrate making and refining that approximation. 

We need a discrete approximation to the difficult function 
sin^ in the range [0, 9q\. Look at the two functions 9 and sin 9 
after dividing by 9; we are taking out the common big part, the 
topic of Chapter 5. The difficult function becomes (sin 9)/ 9. 
The other function, a straight line, is the simple harmonic ap- 
proximation, and is a useful zeroth approximation. But it does 
not produce any change in period as a function of amplitude 
(since the height of the replacement line is independent of 9q). 
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The next approximation does fixes that problem. Use a 
flat line with height (sin0o)/#o- This line is the minimum 
height of (sin0)/9. Why is that choice an improvement on the 
first approximation, using the maximum height of 1? Because 
in this choice, the height varies with amplitude, so the period 
varies with amplitude: This choice explains a physical effect 
that the first choice approximated into oblivion. In this second 
approximation, the torque term (g/ 1) sin 9 becomes 




I 9 ■ 

Starting from the simple-harmonic approximation, this choice is equivalent 
to replacing gravity by a slightly weaker gravity: 



9 -> 9 x 



sin Uq 



9< 



o 



The Taylor series for sin gives 



The fake g is then 



sm e a _ el 

~ 6' 



5fakc = g ( l - ) ■ 



Using this fake g, the period becomes 



2tt* 



fffako 



— 1/2 

To compute g ia ^ c requires another Taylor series: 



(1 + x) 



-1/2 



1 - 



Then 
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This period is an overestimate because it assumed the 
weakest torque adjustment factor: scaling the torque by the 
value of (sin 9) /0 at the endpoints of the swing when 9 = ±# - 
The next approximation comes from using an intermediate 
height for the replacement line. Equivalently, say that the 
pendulum spends half its flight acting like a spring, where the 
torque contains just 9; and half its flight where the torque has 
the term 0(sin#o)/$o- Then the period is an average of the 
simple-harmonic period T — 1~K^/Tjg with the preceding underestimate: 




The next step - and here I am pushing this method per- 
haps farther than is justified - is to notice that the pendulum 
spends most of its time where it moves the slowest. So it 
spends most of time near the endpoints of the swings, where 
the simple-harmonic approximation is the least accurate. So 
the endpoint-only underestimate estimate for T should be weighted 
slightly more than the simple-harmonic overestimate. The 
most recent estimate weighted these pieces equally. To im- 
prove it, count the endpoint estimate, say, twice and the center estimate 
once. This recipe has a further justification in that there are two endpoints 
and only one center! Then the period becomes 




The true coefficient, which comes from doing a power-series solution, is 1/16 
so this final weighted estimate is very accurate! 





3.7 What have you learnt 

Discretization makes hard problems simple. The recipe is to replace a com- 
plicated function by a rectangle. The art is in choosing the height and width 
of the rectangle, and you saw two recipes. In both, the height is the max- 
imum of the original function. In the first, easier recipe, the width is the 
range over which the function changes by a factor of e; this recipe is useful 
for linear exponential decays. The second recipe, the FWHM, is useful for 
messy functions like spectroscope absorption peaks and Gaussians. In that 
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recipe, the width is the width over which the function goes from one-half the 
maximum and then returns to that value. 
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